The spatial propagation and increasing dominance of Gilbertiodendron dewevrei (Fabaceae) in the eastern Congo basin

Though substantial research has been conducted on possible historical, physiological, and symbiotic mechanisms that permit monodominance to occur within tropical lowland rainforests, less is known about the successional rates at which monodominance exerts itself on surrounding forest structures. Here we extend efforts to evaluate the longitudinal dynamics of Gilbertiodendron dewevrei-dominated forest in Central Africa by considering this species’ spatial dynamics. Using three 10-ha censused field plots measured across three time periods, we present the first quantitative estimates of the spatial propagation of Gilbertiodendron into adjacent mixed species forest. Using three analytical strategies, we demonstrate that Gilbertiodendron is increasing in dominance and that monodominant forest patches are expanding into the surrounding forest at a statistically significant rate. The rates of successional advance vary by patch and direction, but average 0.31 m year−1, with speeds greatest in the direction of the prevailing winds. We show that the advancement of Gilbertiodendron is significantly slower than documented rates from other forest ecotones across Central Africa. When paired with stress tolerance traits and ectomycorrhizal associations, these findings help to clarify the means by which Gilbertiodendron dewevrei gains dominance in otherwise species-diverse regions.

Though Gilbertiodendron-dominanted forests have been studied for the better part of a century [2,14,16,20,26], we still lack a unified understanding of the mechanistic processes by which their monodominance occurs, and a clear understanding of their rates of forest succession. To better understand the manner in which a mixed-species forest composition yields to monodominance, Glick et al. [13] examined shifting forest types in the Ituri region of the DRC using a successional trajectory defined by: mixed species forest, in which Gilbertiodendron represented <5% of the total living basal area (BA); transitional forest, in which Gilbertiodendron represented �5% and <50% of the BA; and Gilbertiodendron-dominated forest, in which Gilbertiodendron represented �50% of the BA. Statistically significant results suggest that the forests of the Ituri region have followed a uni-directional successional trajectory leading from mixed species to transitional to Gilbertiodendron-dominated forest. Thomas et al. [27] also presented findings that considered the transition from mixed species forest to Gilbertiodendron-dominated forest. Using two of the three censuses utilized in the present work, these authors found that in the Ituri region of the DRC, the boundary between Gilbertiodendron-dominated forest and mixed species forest was moving towards the latter at a rate of 0.34 m per year (per [28]). Our research was motivated by this initial effort.
In the present study, we sought to extend these efforts by focusing on one over-arching question: Is Gilbertiodendron dewevrei expanding into the surrounding mixed species forest, effectively converting mixed species assemblages into monodominant forest? We focused on three hypotheses that (a) targeted complementary measures of Gilbertiodendron's propagation and displacement over space, along with the associated displacement of mixed species forest; and that (b) capitalized on unique computational strategies. We hypothesized that: 1. When looking across each field site, the spatial extent of Gilbertiodendron-dominated forest patches has increased substantially over time.
2. Relative to patch centroids, Gilbertiodendron-dominated forest patch boundaries have migrated outward, reflecting patch expansion across time. 3. The magnitudes and directions of change over time in Gilbertiodendron BA intensity (BA/ m 2 ), reflect both an increase in BA and greater rates of change around the margins of regions with higher BA concentrations.
Late successional monodominance is generally characterized by slow-growing, shade tolerant species with large seeds and poor dispersal [9]. Slow growth and poor seed dispersal would tend to contraindicate widespread success and short-term influence. However, previous research has indicated that Gilbertiodendron can be a key determinant in forest structure [15,18,21] with increasing success over relatively short time periods [13]. Here, we take a more direct approach to (a) examine whether Gilbertiodendron is advancing into surrounding forest structures, and to (b) quantify the rate of that movement.

Materials and methods
All data management, analysis, and visualization tasks were performed using R v. 3.5.3 (R Core Team, 2019). Specific packages and workflows are described below, and are available in Supporting Information S1 Data.

Field plots and data structures
This research was based on three of four field plots associated with the Ituri Forest Dynamics Project, which have described in detail by [13,29] and references therein. Established in 1994 by the Centre de Formation et de Recherche en Conservation Forestière and jointly managed with the Smithsonian Institution's Center for Tropical Forest Science (CTFS), the Wildlife Conservation Society (WCS), and the Okapi Faunal Reserve (OFR), the Project includes four 10-ha (200 m × 500 m) field plots covering the range of local forest types: mixed species forest, Gilbertiodendron-dominated forest, transitional forest between mixed and monodominant, and forest associated with riparian zones or other microsites. The field plots occur in two groups named for the Edoro and Lenda Rivers: edodo1, edoro2, lenda1, and lenda2. The edoro2 field site was originally created as a control plot. At the time of installation, it was composed entirely of mixed species forest and did not contain any Gilbertiodendron. For this reason, it was omitted from the present analysis. The remaining field plots fall on acidic, welldrained oxisols at 700-850 m above sea level, with limited topographic variability (mean of 18 m between lowest and highest points) [4,29] (Fig 1).
Of particular interest here, the Ituri Forest Dynamics data includes repeated measurements for all stems � 1 cm DBH (diameter at 1.3 m height). From Feb. 2, 1994-Jul. 15 1996, Feb. 2, 2001-Jan. 9, 2002, and Jan. 10-Nov. 30, 2007, three separate censuses were carried out in which every stem was tagged, identified (> 99% to the species level), measured for DBH, and geographically referenced to a local origin. A fourth census was conducted between Aug. 1, 2015 and Apr. 17, 2018, but at the time of this analysis, the data had not been sufficiently cleaned and unified to permit its inclusion in the present work. Our full tabular dataset included 959,312 individual stem measurements from 342,893 unique trees across the four field sites and three censuses. Of these, we relied on 34,499 records of living Gilbertiodendron with valid DBH or BA measurements, and omitted 1,411 records of dead stems lacking suitable diameter measurements. All field plots with Gilbertiodendron contained small diameter seedlings and saplings (� 2 cm DBH) at each census, though the number of small stems was an order of magnitude greater on the lenda field plots (n = 1151-1910 per census), than on the edoro1 plot (n = 158-182 per census). Where Gilbertiodendron stands display synchronous supra-annual masting every 3-5 years [4,18], the small stems might appear to implicate local masting during inter-census years. However, Gilbertiodendron growth is quite slow in shaded understory conditions, and small stems are generally older than they appear [18].
For those of our analysis that required locational information, it was necessary to record geospatial (x,y) locations for each stem. For the 782 records (2%) that lacked coordinate Overview of study area and Ituri Forest Dynamics Project field plots. The location of Ituri field plots in the DRC (a); regional context, showing the two sets of field plots, named for their proximity to the Edoro and Lenda Rivers (b); and local context, showing Edoro (c) and Lenda (d) field plots with minor plot-level hydrological features digitized from [29]. The eastern Edoro plot was omitted from analysis due to lack of Gilbertiodendron. Basemap from Global Land Analysis & Discovery laboratory 2019 Landsat mosaic. Regional hydrological features derived through hydrological analysis of 3 arc second Shuttle Radar Topography Mission digital elevation tile North 01 East 28. https://doi.org/10.1371/journal.pone.0275519.g001

PLOS ONE
The spatial propagation and increasing dominance of a tropical monodominant information due to issues with record keeping, we were able to identify each stem's 20 m × 20 m parent quadrat (a function of the field study design), and then randomly situated each stem within that quadrat. Randomization was performed across both x and y dimensions, and ensured that no two stems were located at the same location. This technique provided unbiased coordinates in aggregate, and suited our needs given the spatial distribution of coordinate-deficient records and the suite of techniques used in subsequent analyses. To move stem records from a local geospatial referencing system to global coordinate system, spatial field plot boundaries were defined as 200 x × 500 m rectangles based on plot origins, which were themselves established using a civilian-grade GPS transceiver. We initially defined the plot corners with respect to the World Geodetic System of 1984 (WGS84), and then projected them to Universal Transverse Mercator Zone 35 North (UTM 35N, m). To covert each stem's original coordinates into the projected coordinate system, we then applied stem-specific offsets relative to each plot origin.

Monodominant patch delineation
To model the translation of Gilbertiodendron-dominated forest over space, we employed three distinct strategies. In the first, we considered changes in the area-weighted mean radius of gyration over both space and time. In the second, we used linear models to consider changes in the projected coordinates of patch boundaries over space and time. In the third, we used vector fields to visualize and quantify the magnitudes and directions of change in Gilbertiodendron basal area (BA) intensity over space and time (Fig 2).
To understand how monodominant forests were changing, we first isolated records from each field site and census that contained Gilbertiodendron dewevrei (n = 9, as edoro2 had no Gilbertiodendron by design). For each of these subsets, we used the spatstat package [30,31] to develop a continuous surface of aggregated Gilbertiodendron BA within 1 m 2 grid cells. Prior to applying a moving-window function to each grid, we implemented an edge correction strategy drawn from the sampling literature. When measuring the contents of a field plot that borders the parent domain, there are implicit edge effects from the portions of plots that extend beyond the measurable domain. Gregoire and Valentine [32] discuss the "reflection method" as one method to account for these effects. In the reflection method, a direction vector from the plot origin to the domain boundary is established, and then "the direction vector is folded back over itself at the boundary and the 'reflected plot' is installed where the vector terminates inside of" the initial plot. Here, we employed the reflection method, but instead of reflecting the plot-equivalent here, to each of the moving windows described below-to make up for the portion that falls outside domain boundary, we reflected the BA intensity surface. Reflection was performed in each of the cardinal and inter-cardinal directions to produce, for each field site, a continuous grid 200 m larger than the original on all sides.
To the expanded BA surface, we used the raster package [33] to apply a circular movingwindow summation with radius 56.42 m, capturing the per hectare BA within the focal neighborhood [34] of each cell. To these grids, a 3 × 3 cell moving window mean was applied to reduce the sinuousity and vertex complexity of the margins of forest patches (described below). These grids, one for each of three sites and two time periods (census 1 and census 3 only), were then cropped to their original dimensions. From these surfaces of BA per hectare, we identified Gilbertiodendron-laden forest patches for use in subsequent analyses. From the BA intensity surfaces, we generated isolines corresponding to 20 m 2 of BA on the lenda1 and lenda2 field sites, and the isolines corresponding to 10 m 2 of BA on the edoro1 field site (since it had substantially less well-developed Gilbertiodendron dominance compared to lenda1 and lenda2). While the Ituri rainforest exemplifies monodominant and mixed species forest assemblages, at times with clear boundaries [16,22], the interspersion and juxtaposition of Gilbertiodendron-dominated forest against other forest types is not so clearly defined within our field data. We selected these 20 m 2 and 10 m 2 thresholds based on an approximate, but ultimately subjective, assessment of what constituted patch boundaries. Finally, using functions from the raster package [33], isolines were converted into discrete (singlepart) polygonal features, adopting the boundary of the field site in cases where the two intersected. These features represented the Gilbertiodendron-dominated forest patches of interest.

Measuring area-weighted mean radius of gyration
Our first approach for quantifying the movement of Gilbertiodendron-dominated forest patches was based on the area-weighted mean radius of gyration [35,36] (1), which we computed for each field site and census. Also known as correlation length, this metric is often used to characterize animal movement. Keitt et al. [35] described it as the mean distance an individual (animal or otherwise) is able to move before reaching a barrier, when placed randomly within suitable habitat. Taking an area-weighted average of radius of gyration values essentially captures the average distance an individual might travel while remaining within the same patch [36]. While trees cannot move, we can treat offspring from a given tree as future versions (a) Raw tabular data was cleaned and coerced into spatial format (i). Basal area was aggregated to grid cells (ii), whose values were minimally smoothed after accounting for edge effects (iii). Monodominant forest patch delineation was based on isolines of BA intensity (iv). Three different different analytical approaches were then used to evaluate spatial propagation of monodominant forest on each field plot (b-d). (b) Scale-independent radius of gyration was computed for each patch (i), and a patch area-weighted mean was computed for each field plot (ii). The final metrics were used to evaluate change between the first and third censuses. (c) Patch vertices from the first and third censuses were paired based on measures of proximity (i). Regression was used to look at changes in latitude and longitude for each of four cardinal quadrants (i). Permutation testing was used to determine if the linking algorithm was producing demonstrable, non-random measures of change (ii). (d) The difference in basal area intensity between the first and third censuses was treated as a topographic surface from which vector magnitude (slope) and direction (aspect) were derived (i). Spatial bootstrapping of the pixelated difference surfaces was used to evaluate plot-level gains or losses across time (ii).
https://doi.org/10.1371/journal.pone.0275519.g002 of the seed-producing parent. This approach provided a speculative way to answer the question, if trees could move, where would we expect tree x to be located n years from now?
Area-weighted mean radius of gyration is directly related to landscape fragmentation, for as values increase, so too does connectivity of suitable habitat. In our usage, "barrier" is synonymous with the edge of a target patch, and by taking an area-weighted mean across each field site, we establish an intuitive measure of the site-wise spatial extent of patches. There are a number of metrics that seek to capture the spatial extent of land cover patches [36], and we elected to use area-weighted mean radius of gyration to avoid scale dependence. For each field site, we subtracted the area-weighted mean values associated with the first census from those of third census in order to understand, for each field site, the average annual landscape-level expansion or recession distance (in meters) of Gilbertiodendron-dominated patches. Radius of gyration is computed as where h ijr is the distance between cell r within patch ij and the centroid of patch ij, using cellcentroid-to-cell-centroid distance, and where z is the number of cells in patch ij. In practice, for each of our six collections of forest patches (three sites × two censuses), we computed the radius of gyration for each conterminous patch before using the patch areas to obtain an areaweighted mean. Patches as small as 1 m 2 in size and containing as few as a single stem or coppice were considered, and an 8-cell ("queen") adjacency rule was employed.

Modeling changes in patch boundaries
In our second strategy for characterizing the movement of Gilbertiodendron-dominated forest patches, we used an algorithm to pair and trace patch boundary vertices over time and space, emulating organism mobility. We used the raster [33] and rgeos [37] packages to isolate only those patches � 300 m 2 , and those whose nearest boundary segments were � 30 m from another patch. Since we were interested in established forest patches, the former rule sought to mitigate the impacts of spurious individuals and small patches of Gilbertiodendron that may materialize and dematerialize over time through natural successional dynamics or disturbance processes. The latter rule was imposed in order to identify unique (i.e., spatially isolated) patches. Species in the Gilbertiodendron genus tend to have deep crowns with poor seed dispersal [23]. The heavy seeds are ballistically dispersed upon dehiscence of the seed pods, limiting advancement across the landscape [18,38]. Hart [18] identified seed dispersal distance to be 5-10 m beyond the vertical projection of the crown edge, and van der Burgt et al. [39] (see also [40]) speculated that the maximum distance "is probably in the range of 50-70 m" (p. 56). Van der Burgt [38] (p. 29) also speculated that the upper bound on dispersal distance from the tallest tree in his study was "likely around 60 m". In light of our limited data, we elected to apply a conservative threshold of 30 m, beyond which we considered target features to be part of a unique patch. Patches that met these criteria were assigned unique identifiers in order to track them over time.
For each unique patch, we employed a matching algorithm to measure the distance between patch boundaries at the first and third censuses. Working with patch boundary vertices (as opposed to inter-vertex line segments, or edges), we matched each vertex from the third census to the nearest vertex from the same patch in the first census in order to link vertices based on the most likely movement trajectory for each point on the patch boundary. Here, "most likely" was equated with the shortest Euclidean distance, which was, in most cases, approximately orthogonal to the boundary curve.
For each unique patch, a first pass of our vertex-matching algorithm moved sequentially through each vertex in the census 3 patch boundary, identifying the closest vertices within the census 1 boundary. The census 3 vertices were then ordered by distance to their matching vertices, and the pairing algorithm was applied to all census 3 vertices in descending order of distance. Because Gilbertiodendron patches generally expanded over time, there were more vertices in the census 3 patch boundaries than in their census 1 counterparts. For this reason, we permitted up to two vertices in the third census to be paired with the same vertex in the first census on the lenda1 and lenda2 field sites (� 2 : 1), and, because of even greater expansion, up to three vertices from third census on the edoro1 field site (� 3 : 1).
For each unique patch in each field site, paired vertices were grouped according to their cardinal quadrants relative to the census 1 patch centroid. Pairings were assigned to east or west quadrants based on a north-south longitudinal meridian cutting through the patch centroid, and to north or south quadrants based on an east-west parallel cutting through the patch centroid. Defining locational quadrants in this way enabled us to use weighted simple linear regression (nlme::gls [41]) to model the change in projected coordinates of the paired vertices as a function of time, with one model for each quadrant, for each target patch, and for each landscape (n = 4 directions × 5 unique patches in total = 20 separate regressions). The regression models we employed can be expressed as, where Y = latitude or longitude, i = 1, . . ., n indexes (pairs of) vertices associated with the third census, j = 1, . . ., 5 indexes unique patches, and q = N, E, S, W captures the directional quadrant. Here, we employ a sub-model for the variance using nlme::varIdent [41,42]. In the generic form, this enables the model to estimate a different σ 2 for each level of a stratification variable s, indexed with s = 1, . . ., S. To paraphrase Pinheiro and Bates [41] (p. 209), this variance models uses S + 1 parameters to represent the S variances, requiring that we impose restrictions on the δ parameter to ensure identifiability. In all cases, δ > 0, δ 1 = 1, and each subsequent δ value (2, . . ., S) represents the ratio of the standard deviations of the sth stratum and the first stratum. In the present application, we have only two strata; s = 1 for census 1 and s = 2 for census 3. After applying (2) to each quadrant of each unique patch, we extracted the model coefficients in order to identify the linear distance of change (meters) in the Gilbertiodendron patch boundary per year.
We employed case resampling-based bootstrapping (n = 10, 000) to compute estimates of mean annual change and the associated bias-corrected, accelerated 95% confidence intervals [43] (boot::boot.ci). Since our pairing algorithm was based more on ecological principles and common logic than common analytical strategies, we also used permutation testing to ascertain whether the results obtained were the product of a non-random and potentially meaningful process. The version of the pairing algorithm we implemented in our permutation tests was identical to that described above, except that we held the vertices from the third census constant and randomly permuted those of the first census. Thus, matched pairs were not necessarily the geographically closest. We computed 10,000 random permutations for each patch within each field site using cluster-based parallel computing. The results were aggregated and the distributions of coefficients for the year predictor were used to evaluate the statistical significance of our initial findings.
Quantifying the flow of basal area. Our third approach for characterizing the movement of Gilbertiodendron-dominated forest patches utilized vector fields. We used spatstat [30,31] to convert our original stem locations to ppp objects (spatially bounded two-dimensional point patterns) with BA as marks, and then subjected these ppp objects to kernel intensity estimation using density.ppp.
We employed a multivariate Gaussian kernel, the Jones boundary correction [44], and a 25 m bandwidth. While we created BA intensity surfaces for our radius of gyration calculations, these relied on the empirical BA ha −1 around each pixel and were therefore insufficiently smooth for vector field analysis. Using the kernel smoothed intensity surfaces, we subtracted that of the first census from that of the third to obtain a difference grid. The difference grid is a three dimensional topographic surface to which we applied slope and aspect calculations for each 1 m 2 cell on the landscape. The slope and aspect dimensions represent vector magnitude and direction, respectively. Using custom R functions, we plotted vector fields of change in BA intensity for each of the three field sites. To enhance interpretative clarity we uniformly scaled the vector magnitudes and visualized a systematic sample of vectors to avoid graphic overcrowding.
In addition to the vector field plots, we utilized a bootstrapping procedure to estimate the mean change in BA ha −1 . To each of the three difference grids, we applied the same reflectionbased edge correction method discussed above [32]. To the expanded grids we applied a circular moving-window summation (radius 56.42 m) in order to capture the empirical sum of the change in BA in 1 ha neighborhoods surrounding each pixel. Each of these grids was then reduced to its original size before using pixel-wise case resampling-based boostrapping (n = 10,000) to obtain a mean and 95% confidence interval for change in BA ha −1 within each field site.

Results
When calculating differences in the area-weighted mean radius of gyration for different field sites containing Gilbertiodendron, we found that all patches were expanding ( Table 1). The edoro1 field site was characterized by younger, less well established patches of Gilbertiodendron that occupied between 11.1% (census 1) and 15.1% (census 3) of the landscape. Over the 12-year study, these patches expanded at a rate of 0.44 m per year. Lenda1 had a single wellestablished patch that covered between 65.6% and 68.0% of the site, a 3.5% change equating to a patch expansion rate of 0.27 m per year. Lenda2 presented two well-defined patches occupying the majority of the site (�82%). Advancement over time was positive but negligible, with a growth rate of 0.03 m per year.
Our isoline movement models generated comparable, if more nuanced, results. Fig 3 visually illustrates changes in Gilbertiodendron-dominated forest patch boundaries over the 12-year study period. Results were most pronounced in edodo1, whose patches were less mature than those of lenda1 or lenda2. The lenda2 field site showed minor recession towards the western margin of the site, but clear expansion into the mixed species forest median that formerly separated two distinct Gilbertiodendron patches in the north. The lenda1 field site showed minor expansion in the northeast, southeast, and southwest, but otherwise exhibited a more steady state than the other two sites. Along the largely coincident patch margins of lenda1 was a natural fluctuation with both expansion and contraction on the order of meters to tens of meters. Table 2 presents the numerical results of our isoline movement modeling. Here, our expectation was that all patches would be expanding at statistically significant rates. Because we were modeling coordinates, and because all of the examined patches fell to the north and east of the UTM 35N coordinate system origin, patch expansion was characterized by increasing latitude in the northern quadrant (positive estimated coefficients, or + in the Expect. column), increasing longitude in the eastern quadrant (+), decreasing latitude in the southern quadrant (-), and decreasing longitude in the western quadrant (-). In aggregate, we found that 88% (7 of 8) of quadrants in edoro1 were expanding in accordance with expectations. This number is reduced to 75% (3 of 4) for lenda1 and to 63% (5 of 8) for lenda2. Thus, 75% (15 of 20) of quadrants reflect aggregate expansion of Gilbertiodendron into areas of (previously) lesser Gilbertiodendron BA intensity. The lenda1 field site did not exhibit any statistically significant changes, as can be seen in Fig 3. Meanwhile, 75% (6 of 8) of quadrants in edoro1 and 50% (4 of 8) of quadrants in lenda2 presented statistically significant patch boundary movement. Of the changes that were significant, only two (edoro1 patch 1 N, lenda2 patch 1 E) reflected recession instead of expansion.
Since the prevailing wind patterns in the Ituri province are from the northeast, we might hypothesize that we would see greater expansion in the southerly and westerly directions than we would in northerly or easterly directions if the trajectories of ballistically dispersed seeds were being influenced by the wind. At both edoro1 patch 1 and lenda2 patch 2, the greatest rates of expansion were in the south and west quadrants. In edoro1 patch 2, the second and third greatest rates were in the south and west quadrants, and in lenda2 patch 1, the second greatest rate was in the south quadrant (the west had negligible change). The sole patch at lenda1 showed mixed results. When averaging across like quadrants from different field sites, the greatest rates of expansion were to the south and west (Table 2).
With few exceptions, our permutation tests (S1-S3 Figs) demonstrate that the algorithm we used to associate Gilbertiodendron-dominated forest patch boundaries in the first census to those in the third, identified non-random and potentially meaningful associations. The only test for which our original modeled results were not highly significant (p � 0.0008) was that for the western quadrant of lenda2 patch 2. The probability of obtaining this particular model's estimated coefficient was 0.06. This patch was the smaller of the two lenda2 patches, and the limited number of vertices present in the isoline boundary may have led to this marginally significant result. Fig 4 presents the results of our vector field analysis, illustrating changes in the BA for each 1 m 2 pixel between the first and third censuses. As expected, the vast majority of changes are positive, reflecting a net increase in BA over time. This is empirically captured in our bootstrap results: the edoro1 field site increased by an average of 0.35 ± 0.002 m 2 BA ha −1 , the lenda1 field site increased by an average of 0.84 ± 0.004 m 2 BA ha −1 , and the lenda2 field site increased by an average of 0.41 ±0.01 m 2 BA ha −1 . Areas of change contain larger magnitudes along the peripheries of high and low concentrations regions. This is reflected in shorter length vectors over the peaks and valleys in Fig 4, with longer length vectors on the slopes leading to and away from the valleys and peaks. We would expect the greatest rates of change to occur just beyond the isoline derived from the third census' data, as these are the areas in most rapid transition. We would also expect vector directions to reflect an intuitive expansion away from the heart of each patch. The edoro1 field site provides the clearest example of this behavior (Fig 4a lower panel). The southern margin of the sole patch at lenda1 contains a band of increasing BA that parallels the expanding isolines (Fig 4b lower panel), and the northernmost patch at lenda2 also reflects positive gains and expansion along the patch margin (Fig 4c  lower panel).

PLOS ONE
The spatial propagation and increasing dominance of a tropical monodominant The detailed plots in the lower portion of Fig 4 reflect moderate agreement between the expanding isolines emphasized in Fig 3 and the directions and magnitudes of change. Discrepancies between these two analytical approaches can be traced to the underlying data. Our isoline results are based on minimally smoothed surfaces of empirical BA/ha (i.e., a uniform kernel), while our vector field results are based on kernel smoothed BA intensity surfaces built with Gaussian kernels and a 25 m kernel bandwidth. It is well known that kernel smoothing is sensitive to changes in the bandwidth parameter [30,45], and we caution the reader in drawing too close an association between our isoline and vector field results. Kernel intensity estimation is generally regarded as a non-parametric technique for exploratory data analysis, and we offer our vector field findings primarily as useful heuristics.

Discussion
Our three sets of results are in concordance and collectively describe a successional process in which Gilbertiodendron dewevrei infiltrates and displaces adjacent closed-canopy mixed species forest, ultimately leading to the spatial propagation and increasing presence of monodominant forest. Of particular interest here is the rate at which Gilbertiodendron is expanding into mixed species forest. To date, there has been limited research on the successional rates of closed-canopy tropical forests with monodominant climax phases, making it challenging to place our findings into a larger context. As points of comparison, we first draw on a number of studies examining tropical monodominant-mixed species interactions that still exist ( § Monodominant-mixed species interactions), before considering successional modes and advancement rates derived from other ecotones ( § Successional modes and rates of forest advance). We then consider gap dynamics in patch interiors, as they relate to Gilbertiodendron's spatial propagation ( § Invasibility and gap dynamics), and end by highlighting the inherent limitations in our findings ( § Limitations).

Monodominant-mixed species interactions
The idea that monodominant forest species possess unique characteristics that enable them to invade adjacent mixed species forests, arose nearly a century ago, if not earlier [46]. The Mora excelsa-dominated forests of Trinidad, which continue to invade the nearby Crappo-Guatecare forest, were an early focus [2,[46][47][48][49]. Mora excelsa's high canopy presence [2], relatively sharp monodominant-mixed species boundaries [2], and dense sapling growth within colonization fronts [46], share parallels with the forests at Ituri. Eggeling [50] also presents an early work, in which he describes the slow invasion of a mixed species forest by Cynometra alexandri in the Budongo rainforest of Uganda. As at Ituri, the author describes several forest types that represent a successional trajectory, placing the mixed species forest as seral to the Cynometra-dominated climax forest (c.f., [13]). Eggeling [50] also describes how Cynometra appears before the mixed species forest has peaked, increasing pressure on other species and slowly asserting itself through understory growth, not unlike what has been observed with Gilbertiodendron [17,21]. From these early works, Beard [2] provides us with a rare map (c.f. Fig 6) of speculated Mora excelsa advancement trajectories into mixed species forest.
Hart et al. [16] write of the same Ituri rainforest examined in the present study, describing the association between Gilbertiodendron-dominated forest and mixed species forest in relation to other tropical monodominant forest complexes. In line with Eggleling's [50] observations, Hart et al. [16] present a case for mixed species forest as seral to Gilbertiodendrondominated forest, and they view monodominant forest boundaries as a reflection of the slow dispersal of the monodominant species following reduction in its historical range by largescale disturbance. They did not find any reason to think Gilbertiodendron, once established, would not continue to assert its dominance over multiple generations (i.e., Connell and Lowman's [8] type I monodominance). Makana et al. [21] present similar findings for the Ituri sites.
Henkel [7] describes the life history and growth traits of Dicymbe corymbosa, a Caesalpiniaceae monodominant in the mountains of western Guyana. He identified a number of structural traits that mirror those at Ituri and which suggest a similar, though perhaps slower, set of successional dynamics. These traits include reduced stem densities, reduced species diversity, and increased basal area in monodominant plots relative to mixed species plots, as well as the complete absence of Dicymbe in mixed species plots. Like Gilbertiodendron, the growth characteristics of Dicymbe accentuate a type I monodominant [8] with a strong ability to gain and retain dominance through mast fruiting, high recruitment rates, litter trapping, an ectomyccorhzial association, conspecific suppression through shade, and a coppicing growth form that leads to "cumulative reiteration" [7,51]. In nearby French Guiana, Charles-Dominique et al. [52] present a study on the colonization front of an understory palm, Astrocaryum sciophilum, at its interface with undisturbed lowland rainforest. Simulation studies revealed a clear pattern of heterogeneous advance, and field observations were used to estimate an advancement rate of 2.3 m year −1 for Astrocaryum into mixed species forest. To our knowledge, this is the only published rate of advancement of a colonizing tree (or tree-like) species within closed-canopy tropical rainforest.
Marshall [46,47], Gunther [48], Beard [2], Eggeling [50], Hart et al. [16], and Glick et al. [13] all present examples in which mixed species forest is seral to monodominant forest. In contrast, Peh et al. [23] test the notion that Gilbertiodendron monodominance is seral to mixed species forest, perhaps building from Hart et al.'s [16] discussion of theoretical patterns of succession among adjacent monodominant and mixed species forest communities. Although they identified seven mixed species forest species succeeding in Gilbertiodendron's understory, this is a relatively small percentage in relation to the 438 unique species at Ituri [13]. While we agree with Peh et al.'s [23] assertion that some species can establish themselves in Gilbertiodendron-dominated forest, we speculate that this may be a function of the transitional forest structure that exists at the interface between mixed species forest and monodominant forest [13]. Further, if Gilbertiodendron-dominated forest was a successional community, we would generally expect the species to possess the characteristics of an early successional pioneer [16]. Yet, to paraphrase Hart et al. [16], we have little to no evidence of a more-specialized species that can or will succeed Gilbertiodendron in the absence of major disturbance.
Similarities between different monodominant forests have been examined previously [16], and the above-referenced studies indicate that monodominant-mixed species associations also share a suite of characteristics. As at Ituri, transitional forest structures, limited monodominant presence in mixed species forest, climax monodominance, colonization fronts, shade tolerance, and the ability to retain dominance over more than one generation, are characteristics present in comparable forest associations. Only Peh et al. [23] provide any evidence of successional processes that question Gilbertiodendron-dominated forest as a climax state. Since only one of the above-referenced studies provides us with an estimated rate of a species' advance into mixed species forest (neither a dominant nor canopy species), we now turn to forest/nonforest ecotones for additional insight.

Successional modes and rates of forest advance
Although forest/non-forest ecotones share markedly different successional trajectories than those of the closed-canopy forest at Ituri, they provide some of the very few estimates of forest encroachment rates that we could locate. Forest/non-forest ecotones follow two primary modes of forest succession: (a) linear edge diffusion [53], in which succession occurs as an advancement front oriented roughly perpendicular to the progression trajectory; and (b) a more stochastic process of clump formation surrounding a nucleus of opportunistically established early successional flora, followed in turn by the coalescence of clumps that gain adjacency through edge diffusion [54,55].
References to clump coalescent succession first appeared in the early 1960s, if not earlier. Hopkins [56] conjectured that rapidly spreading forest overtaking savanna in western Nigeria was a function of isolated forest foci all spreading out into the surrounding landscape. By comparing two hand-drawn maps (with error acknowledged by the author), he estimated the exceptional annual expansion rate of 33 m year −1 . Whittaker et al. [57] and Archer [54, 58] described succession via opportunistic clump formation as a two-phase process in which there was an open continuous phase characterized by a grassy matrix, and a discontinuous phase in which the matrix was punctuated by clumps of woody vegetation. They examined a mesic grassland/savanna/forest ecotone in southern Texas, USA. (c.f. [58] Fig 1) for which the mean absolute growth rate across clusters of all sizes ranged from 1.8 to 3.2 m year −1 radially [54]. In their carbon isotope-based analysis of a savanna/forest ecotone in east-central Cameroon, Guillet et al.
[59] described clump succession as non-linear and faster than edge diffusion. They felt that clump coalescence was the dominant successional process in their study area, and that it was historically responsible for relatively rapid colonization of the savanna. Puyravaud et al. [60] observed that vegetative thickets served as nodal centers in a successional process leading from grassland to rainforest in the Western Ghats of southern India and, if left undisturbed, thickets would coalesce to form secondary rainforest. Favier et al.
[53] studied forest/savanna ecotones in coastal Republic of Congo and found both modes of succession at play in the conversion from savanna to forest.
Relative to clump succession, succession by edge diffusion may be easier to recognize in the field, and has been discussed many times over the years. Hopkins [56] described the fluctuating boundary of forest and savanna in western Nigeria, noting a forest encroachment rate of 3-4 m year −1 . Citing Miège's [61] work on savanna/forest boundaries in Côte d'Ivoire, Hopkins [62] refers to an expansion rate of 2.5-3.5 m year −1 over a decade or more. Schwartz et al. [63] focused on the Mayombe region of coastal Republic of Congo, and presented a case for at least four successive phases in a linear expansion that converted savanna to sparse secondary forest. In a second study of the same area [64], analysis of vegetation composition, vegetation structure, and soil organic carbon isotope ratios revealed that enclosed savanna was being overtaken by forest at a rate of 0.14 to 0.75 m year −1 , with a mean rate of 0.2 to 0.5 m year −1 . Delègue et al. [65] also used carbon isotope analysis in their work on savanna/forest mosaic of coastal Gabon. They estimated the current rate of encroachment to be approximately 1 m year −1 . Fabing [66] (as cited in [53]) used aerial imagery to estimate that coastal forest was advancing into savanna at a rate of 1 m year −1 in the Republic of Congo. Favier et al.
[53] also studied littoral Republic of Congo and observed that, in moving from forest towards savanna, there were progressively smaller trees and less shade-tolerant species, indicative of a colonization front. They estimated the edge diffusion rate of advancement to be approximately 1 to 2 m year −1 based on a model of diametrical growth of Aucoumea klaineana. Youta Happi [67] reported on the diffusion of forest patches into savanna in Cameroon between 1951 and 1996 (c.f., Fig 48). He found forest to be advancing between 0.6 and 2 m year −1 , with certain forest edges advancing less than 0.25 m year −1 . The author acknowledged the limitations in accuracy, and cautioned the reader in holding too strict an interpretation of these rates. Cuni-Sanchez et al. [68] studied savanna/forest transition and boundary dynamics in Central Gabon, where monodominant Okoume forest was a seral phase following a colonizing forest front. They did not quantify linear rates of advancement of any forest type, but found succession to be slow; after 20 years, not a single field plot could be characterized as the subsequent stage in their putative successional scheme.
We can see from Table 3 that our findings from Ituri are in line with some of the published rates of forest succession in Central Africa. In comparing all of our quadrant models (n = 20) to the range of forest advancement rates from across Central Africa, however, we find that the rates of expansion at Ituri are significantly slower than other sites (-0.75 ± 0.51/-1 m yr −1 , p = 0.01). In particular, there is a much lower expansion rate associated with the eastern Ituri quadrants (0.09 m yr −1 ). This is largely a function of a single outlying value which, when removed, increases the estimated expansion rate and narrows the confidence interval (Table 3). Still, the rates of expansion are significantly slower than other Central African sites (-0.67 ± 0.5/-1 m yr −1 , p = 0.02). Findings from Thomas et al. [27], whose field data is encapsulated by the presented study, have been included in Table 3 for reference. While Thomas et al, [27] looked over a more limited time period and their analytical approach differed from any of those used here, they presented a comparable rate of advance.
Overall, these findings are not unexpected given that, with the exception of Charles-Dominique et al.'s [52] work, the figures in Table 3 were not derived from undisturbed closed-canopy lowland rainforest communities. Our findings agree with previous research in the Ituri region suggesting that the dispersal of Gilbertiodendron is likely to be slow, given this species' slow growth rates and the fact that it has the heaviest seeds of any local species [12,15]. The primary mode of forest succession taking place in our field sites is edge diffusion, as is particularly visible in edoro1 (Fig 3a). However, the northern portion of lenda2 (Fig 3c) implies a clump coalescent successional process may also be taking place. Here, we observe a central nucleus of Gilbertiodendron that is gaining in BA and, in conjunction with the edge diffusion of the primary patches, may eventually lead to a Gilbertiodendron isthmus cutting through the mixed species forest. This type of forest structure, in which islands of one forest type exist within another, was observed by Hart [16] many years ago. Hart [18] observed that the highest rates of Gilbertiodendron seed survival occur just beyond conspecific mast zones, perhaps as a function of increased survival benefits through reduced predator presence. This implies that transitional zones at the intersection of mixed and monodominant forest types are where we are likely to see the greatest expansion of monodominant forest, and that edge diffusion may be a more common mode of succession than clump formation or coalescence. This seems evident from Fig 3 and from the shifting forest classes presented in Glick et al. [13]. Gandhi et al. [55,69] provide an alternative, theoretical lens through which to consider our observations of monodominant patch expansion. They observed a "critical nucleus size" as a spatial size threshold above which a patch of the competitively superior species will continue to grow, and below which it is not apt to persist. They found that the curvature of the interface between interacting species helps dictate the future fate of the species (c.f. [55] Fig 3). Patches tended to take on circular forms over time, and perturbations in the propagation front tended to decay, returning the patch to its more stable circular form. This aligns with their findings on cluster size, as once a cluster reaches the critical nucleus size, its mean perimeter curvature is reduced, placing the nucleating species in the competitively advantageous position within the interaction zone. The larger the patch, the shallower the mean curvature, and the greater the competitive advantage of the superior species. The dynamics described by Gandhi et al. [55] seem at play in our data, and where Gilbertiodendron has a number of advantageous life history traits [17], it appears that the critical nucleus size for the species is quite small (e.g., Fig 3a). The forest structures at Ituri provide a promising location for future field tests of Gandhi et al.'s [55,69] findings.

Invasibility and gap dynamics
The introduction, success, and range expansion of exotic species in a vegetative community is often associated with disturbance, whether natural or human-mediated [70]. This is true even of tropical rainforest communities where canopy gaps are the most common opportunistic channel for changes in species composition. Gap dynamics theory suggests that microclimatic conditions created by canopy openings are critical determinants of succession and invasibility within otherwise contiguous vegetative blocks [71,72]. Yet, a number of studies [2, 23, 46-48, 50, 52], including the present one, suggest that exotic species can establish themselves in mature, closed-canopy tropical (and temperate [73]) forests in the absence of disturbance [23]. This is of interest in that we see both the propagation of Gilbertiodendron into closed canopy mixed species forest (Fig 4), and the formation of gaps in otherwise contiguous forested blocks (Fig 5).
The topographic depressions shown in Fig 4 reflect a general loss in BA across time, and highlight the way in which these landscapes are dynamic and vacillate in cumulative Gilbertiodendron presence. Using an example from the lenda2 field site, Fig 5 seeks to illustrate how gap dynamics actually manifest at Ituri. Here, we note the formation of one clear gap between the first and second census, as well as the formation of another set of gaps between the second and third censuses. Panel (c) of Fig 5 helps to illustrate the gap formation and closure process. The notable decline in BA between the second and third censuses, within the most pronounced gap, is primarily attributable to the loss of two stems (labeled), which were the two largest of those in the focal area. Pronounced changes in plot-level biometrics from select stem loss have been observed in other successional research [68]. In concert with stem loss, panel (c) shows high rates of recruitment (red points) in close proximity to these two lost stems.
There are a number of physiological and life history traits that enable Gilbertiodendron to advance into both closed and open canopy structures. These have been considered in greater

PLOS ONE
The spatial propagation and increasing dominance of a tropical monodominant detail elsewhere [9,13,[15][16][17][18], but it is worth highlighting Gilbertiodendron dewevrei's unique bimodal physiological response to light. Prolific reproduction and understory regeneration speaks to Gilbertiodendron's ability to sustain itself under the deep shade of a closed canopy [18,21], while the species has also shown an unexpectedly large growth response to strong light [19]. So, while the mature forests at Ituri demonstrate the gap dynamics we would expect of mature tropical forests, Gilbertiodendron is unusual in that it may be able to out-compete other late successional species during both the invasion of closed canopy forest and the recolonization of gaps [19].

Limitations
One of the fundamental challenges with examining closed-canopy forest succession is that the forest margins of interest are not defined by clear boundaries such as those in a savanna/forest system. Thus, a choice has to be made in how the forest margin is located, and we recognize that our selection of isolines of 10 m 2 and 20 m 2 of BA introduces a level of abstraction. However, where BA intensity is a continuous spatial process, we argue that comparison of any BA intensity level that does not fall at the extremes of the field site (e.g., not in the lower or upper 10% of BA intensity values) will provide meaningful results. The relationship between tree age and BA accumulation is non-linear and quasi-asymptotic. Depending on what BA intensity is selected, the advancement rates might therefore appear different. However, changing the intensity of interest is equivalent to reframing the focus on the BA intensity associated with a particular stem age class.
Another challenge we face in quantifying the spatial propagation of shifting boundaries revolves around edge effects. Edge effects are well known and we have addressed them here using the reflection method [32] and the Jones boundary correction [44]. However, in our analysis of shifting isolines, we are still faced with the fundamental limits of the plot borders. From an analytical perspective, a Gilbertiodendron-dominated forest patch cannot expand beyond the plot border, effectively biasing our model-based estimates of advancement rates (e.g., Table 2, lenda2 patch 1 E). The solution to this problem would be to establish a longterm monitoring plot in which the spatial extent of the plot captured the entirety of a mixed species forest patch, such that estimates of encroachment of monodominant forest could be quantified using data contained entirely within the spatial domain of the study area. To some extent, our vector field plots (Fig 4) help to overcome this hurdle. They illustrate that, even at the margins of our field sites, BA intensity is increasing. This indirectly reflects the continued propagation of Gilbertiodendron outside of the measured spatial domain.
Finally, the present findings apply to long-term monitoring plots in a particular portion of the DRC. To better understand where our findings fall along a spectrum of common to exceptional Gilbertiodendron movement patterns, similar analyses would need to be conducted with data from other portions of Gilbertiodendron's range. These efforts, drawing on local-level data, would be well supported by the spatially expansive view offered by remote sensing. Given the high cadence and historical span of available data products, analysis of moderate or high resolution multispectral imagery could offer rich insights into the propagation or contraction of larger monodominant patches. As additional high resolution data products enter the open market, we see this as a fruitful area of future research.

Conclusions
The analyses presented here provide direct evidence of the slow advance of Gilbertiodendron dewevrei into mixed species forest in the Ituri region of Central Africa, and provide the first quantitative estimates of the rates of advance of a monodominant species within a closed-canopy tropical forest community. All of the analytical strategies we explored reflect that Gilbertiodendron-dominated forest patches are expanding in area and dominance (increasing BA). Our estimates of Gilbertiodendron's encroachment rate on mixed species forest vary by cardinal direction, with a mean expansion rate of 0.31 m year −1 and greater rates of expansion in the direction of the prevailing winds. Mature patches of Gilbertiodendron present spatial structures typical of mature forest, with gap formation through loss of large stems and regeneration of new Gilbertiodendron recruits. The particular set of spatial analytical strategies we present provides a unique and complementary view of Gilbertiodendron dynamics that is not easily captured through alternative approaches, and we see continued utility in exploring tropical monodominance through a spatial lens. Graphical results of permutation tests (n = 10, 000 iterations) of isoline movement models for site lenda2, patch 1 (a) and patch 2 (b). Histograms capture the null distribution of movement coefficients (m y −1 ) for each directional quadrant: N, S, W, E (clockwise from upper left). The black vertical line is the mean of the distribution, and the red vertical line captures the results of our observed algorithmic pairing of isoline vertices between the third and first censuses. All results highly significant at α=0.05, except patch 2 W (p = 0.06). (TIF) S1 Data. R code for geospatial and statistical analysis. A compressed directory of R data science scripts, custom functions, and support files used to conduct the geospatial and statistical analyses presented herein. (ZIP) Hart, and from Corneille E. N. Ewango. The authors are grateful to Dr. Chad Oliver for his comments on an early version of the manuscript; to Corneille Ewango for his contributions to field work and botanical identification of tree species; to Jacob Madidi, Floribert Bujo, Jacques Mukinzi, Victor Lubuku, and Dieudonné Batido for their assistance with field work; and to Suzanne Lao and Marie-Noelle Ndiu for their assistance with data entry and management.